#include "CalculationKit.h"
#include <iostream>
#include <numeric>
using namespace std;

float Calculationkit::calculateAverage(const vector<float>& data) {
    if (data.empty()) {
        return 0.0f; 
    }

    float sum = accumulate(data.begin(), data.end(), 0.0f);
    return sum / data.size();
}

float Calculationkit::dot_product(float a[3], float b[3]) 
{
    float result = 0;
    for (int i = 0; i < 3; i++) {
        result += a[i] * b[i];
    }
    return result;
}

float Calculationkit::module_length(float a[3]) 
{
    float result = 0;
    for (int i=0; i < 3; i++)
    {
        result += pow(a[i], 2);
    }
    return sqrt(result);
}

void Calculationkit::normalize_vector(float a[3])
{
    float mag_a = module_length(a);
    for (int i=0; i < 3; i++)
    {
        a[i] = a[i]/mag_a;
    }
}

void Calculationkit::mul_vector(float a[3], float length)
{
    for (int i=0; i < 3; i++)
    {
        a[i] = a[i] * length;
    }
}

float Calculationkit::angle_between_vectors(float a[3], float b[3]) 
{
    float dot = dot_product(a, b);
    float mag_a = module_length(a);
    float mag_b = module_length(b);
    float cos_theta = dot / (mag_a * mag_b);
    float theta = acos(cos_theta) * 180 / M_PI;
    // cout << " theta_r: " << theta << "  ";
    return theta;
}

// TRandom3 fRandom;
// float Calculationkit::smearing(float energy_actual, float energyRef, float resolution_at_Eref)
// {   
//     if (energy_actual > 0)
//     {
//         float eRes = resolution_at_Eref * sqrt(energyRef / energy_actual) / 2.35 / 100.0;
//         float gain = fRandom.Gaus(1.0, eRes);
//         float energy_smearing = energy_actual *gain;
//         return energy_smearing;
//     }
//     return 0;
// }


float Calculationkit::camp_calculate(float energy_elec, float energy_scatter)
{   
    energy_elec = energy_elec*1e3;
    energy_scatter = energy_scatter*1e3;
    float cos_em = 1 + Me_MUL_C2/(energy_elec + energy_scatter) - Me_MUL_C2/energy_scatter;
    // cout << " " << Me_MUL_C2/(energy_elec + energy_scatter) << " " << Me_MUL_C2/energy_scatter << " ";
    cos_em = (cos_em > 1.0) ? 1.0 : ((cos_em < -1.0) ? -1.0 : cos_em);
    // cout << "cos_em: " << cos_em << " " << endl;
    float theta = acos(cos_em) * 180 / M_PI;
    // cout << "theta_em: " << theta << "  ";
    return theta;
}


float Calculationkit::calculate_angle_error(float e_vector[3], float electron_energy, float gamma_vector[3], float gamma_energy, float vector_actual[3])
{
    // cout << "electron_energy: " << electron_energy ;
    electron_energy = sqrt(pow(electron_energy, 2) + electron_energy*1.022);
    // cout << " electron_p: " << electron_energy << endl;
    normalize_vector(e_vector);
    cout << "e_vector: (" << e_vector[0] << ", " << e_vector[1] << ", "<< e_vector[2] << ")" << endl;
    mul_vector(e_vector, electron_energy);
    normalize_vector(gamma_vector);
    mul_vector(gamma_vector, gamma_energy);
    float vector_incidence_calculate[3];
    for (int i=0; i < 3; i++)
    {
        vector_incidence_calculate[i] = round((e_vector[i] + gamma_vector[i])*100)/100;
    }
    normalize_vector(vector_incidence_calculate);
    // float angle_error = round(angle_between_vectors(vector_incidence_calculate, vector_actual)*10)/10;
    // if (vector_incidence_calculate[0] < 0) angle_error = angle_error*(-1);
    cout << "vector_incidence_calculate:("  << vector_incidence_calculate[0] << ", "
    << vector_incidence_calculate[1] << ", " << vector_incidence_calculate[2] << ")" << endl;
    // << "angle_error: " << angle_error << endl;
   
    return (vector_incidence_calculate[0] - vector_actual[0]);
    // return angle_error;
}

float Calculationkit::calculate_angle_error_y(float e_vector[3], float electron_energy, float gamma_vector[3], float gamma_energy, float vector_actual[3])
{
    // cout << "electron_energy: " << electron_energy ;
    electron_energy = sqrt(pow(electron_energy, 2) + electron_energy*1.022);
    // cout << " electron_p: " << electron_energy << endl;
    normalize_vector(e_vector);
    cout << "e_vector: (" << e_vector[0] << ", " << e_vector[1] << ", "<< e_vector[2] << ")" << endl;
    mul_vector(e_vector, electron_energy);
    normalize_vector(gamma_vector);
    mul_vector(gamma_vector, gamma_energy);
    float vector_incidence_calculate[3];
    for (int i=0; i < 3; i++)
    {
        vector_incidence_calculate[i] = round((e_vector[i] + gamma_vector[i])*100)/100;
    }
    normalize_vector(vector_incidence_calculate);
    // float angle_error = round(angle_between_vectors(vector_incidence_calculate, vector_actual)*10)/10;
    // if (vector_incidence_calculate[0] < 0) angle_error = angle_error*(-1);
    // cout << "vector_incidence_calculate:("  << vector_incidence_calculate[0] << ", "
    // << vector_incidence_calculate[1] << ", " << vector_incidence_calculate[2] << ")" << endl;
    // << "angle_error: " << angle_error << endl;
   
    return (vector_incidence_calculate[1] - vector_actual[1]);
    // return angle_error;
}


float Calculationkit::calculate_distance_two_hits(float hit0[3], float hit1[3])
{
    float pow_distance = 0.;
    for(int i=0; i<3; i++)
    {
        pow_distance += pow((hit0[i]-hit1[i]), 2);
    }
    return sqrt(pow_distance);
}
